/**
 * @file   potentials.c
 * @author xlhu <xlhu@Nano>
 * @date   Mon May 27 17:32:24 2013
 * 
 * @brief  
 * 
 * 
 */

#include "potentials.h"

/**
 * 1. The following four functions are
 * vector and scalar potentials 
 * that comes from the literatures
 */

void f_xcL0(double complex *fxcL0, double *rho_0, int Nx, int Ny, int Nz)
{
    /******************************************
  longitudinal kernel  f_xcL0 at limit omega=0

  Wigner_Seitz radius: r_s=(3/rho_0/4/pi)^(1/3)
  *******************************************/ 

  int ix=0, iy=0, iz=0;

  double complex fxcl0=0;

  /*************************
  Parameter in VWN
  **************************/
  double A=0, x_0=0, b=0, c=0, x=0, X=0, Q=0, X_0=0, epsilon_x=0, epsilon_c=0, de_x_dr=0, dx_dr=0, dX=0, de_c_dr=0, ddX=0, dde_x_ddr=0, dde_c_ddr=0;
  double de_c_dx=0, ddx_ddr=0, dde_c_ddx=0;;
  //double rho_02=0;

  int icase=0;
  /*************************************
  icase=1, limit \omega \to 0; LDA  VWN
  **************************************/
  icase=1;

  if(icase==1){
    A=0.0310907; x_0=-0.10498; b=3.72744; c=12.9352; Q=sqrt(4.0*c-b*b); X_0=x_0*x_0+b*x_0+c; ddX=2.0;

    for(ix=0;ix<Nx;ix++){
       for(iy=0;iy<Ny;iy++){
          for(iz=0;iz<Nz;iz++){
             //rho_02=rho_0[iz+Nz*(ix+Nx*iy)]*rho_0[iz+Nz*(ix+Nx*iy)]; 
             
             x=pow(3.0/4.0/PI/rho_0[iz+Nz*(ix+Nx*iy)], 1.0/6.0);

             X=x*x+b*x+c;
             
             de_x_dr=-1.0/4.0/PI*pow(3.0*PI*PI,1.0/3.0)/pow(rho_0[iz+Nz*(ix+Nx*iy)], 2.0/3.0);

             dx_dr=-1.0/6.0*pow(3.0/4.0/PI,1.0/6.0)/pow(rho_0[iz+Nz*(ix+Nx*iy)],7.0/6.0);

             dX=2.0*x+b;

             de_c_dx=(A*((2.0*X-x*dX)/(x*X)-4.0*b/((2.0*x+b)*(2.0*x+b)+Q*Q)-b*x_0/X_0 *((2.0*X-(x-x_0)*dX)/((x-x_0)*X)-4.0*(b+2.0*x_0)/((2.0*x+b)*(2.0*x+b)+Q*Q))));
             de_c_dr=de_c_dx*dx_dr;
  
             ddx_ddr=7.0/36.0*pow(3.0/4.0/PI,1.0/6.0)/pow(rho_0[iz+Nz*(ix+Nx*iy)],13.0/6.0);

             dde_x_ddr=1.0/6.0/PI*pow(3.0*PI*PI,1.0/3.0)/pow(rho_0[iz+Nz*(ix+Nx*iy)],5.0/3.0); 

             dde_c_ddx=(A*((dX-x*ddX)/(x*X)-(2.0*X-x*dX)*(X+x*dX)/(x*x*X*X)+16.0*b*(2.0*x+b)/(((2.0*x+b)*(2.0*x+b)+Q*Q)*((2.0*x+b)*(2.0*x+b)+Q*Q))-b*x_0/X_0*((dX-(x-x_0)*ddX)/((x-x_0)*X)-(2.0*X-(x-x_0)*dX)*(X+(x-x_0)*dX)/((x-x_0)*(x-x_0)*X*X)+16.0*(2.0*x_0+b)*(2.0*x+b)/(((2.0*x+b)*(2.0*x+b)+Q*Q)*((2.0*x+b)*(2.0*x+b)+Q*Q)))));

             dde_c_ddr=dde_c_ddx*dx_dr*dx_dr+de_c_dx*ddx_ddr;

             fxcl0=2.0*(de_x_dr+de_c_dr)+ rho_0[iz+Nz*(ix+Nx*iy)]*(dde_x_ddr+dde_c_ddr);

             fxcL0[iz+Nz*(ix+Nx*iy)]=fxcl0; 
          }
       }
    }
  }
  else{
    printf("Wrong reference in eta_xc_xi_xc\n"); getchar();
  }
}

/**
 * longitudinal and transverse kernel  f_xcL,T
 * Wigner_Seitz radius: r_s=(3/rho_0/4/pi)^(1/3)
 * iVK
 *    1: NCT with limit omega \to 0
 *    2: NCT with analytical function
 *    3: QV with analytical function
 **/
void f_xcLT(double complex *f_xcL, double complex *f_xcT, double *rho_0, 
	    double complex omega_input, int iVK, int Nx, int Ny, int Nz)
{
  double omega=0;
  omega=creal(omega_input);

  int icase=iVK;

  FILE  *fout1;
  int ix=0, iy=0, iz=0;
  int ii=0;
  
  double  r_s=0;
  double  x0=0.5, x1=1, x2=2, x3=3, x4=4, x5=5;
  double  *n=(double *)malloc(6*sizeof(double ));

  /*******************
  variables in NCT 0
  *******************/
  double  *mu=(double *)malloc(6*sizeof(double ));
  double  *f_T=(double *)malloc(6*sizeof(double ));

  /*******************
  variables in NCT
  *******************/
  double  *k_F=(double *)malloc(6*sizeof(double));
  double  *epsilon_F=(double *)malloc(6*sizeof(double));
  double  *omega_pl=(double *)malloc(6*sizeof(double));
  double  *fxcL_0=(double *)malloc(6*sizeof(double));
  double  *fxcL_inf=(double *)malloc(6*sizeof(double));
  double  *fxcT_inf=(double *)malloc(6*sizeof(double));
  double  *beta=(double *)malloc(6*sizeof(double));
  double  *c0=(double *)malloc(6*sizeof(double));
  double  *c1=(double *)malloc(6*sizeof(double));
  double  *omega1=(double *)malloc(6*sizeof(double));
  double  *omega2=(double *)malloc(6*sizeof(double));
  double  *d0=(double *)malloc(6*sizeof(double));
  double  *d1=(double *)malloc(6*sizeof(double));
  double omega0=0;
  double g_x=0;
  double  *imag_fxcL=(double *)malloc(6*sizeof(double));
  double  *imag_fxcT=(double *)malloc(6*sizeof(double));
  double  *real_fxcL=(double *)malloc(6*sizeof(double));
  double  *real_fxcT=(double *)malloc(6*sizeof(double));
  double inter=0;
  double inter_imag_fxcL=0;
  double omegaP=0;
  double omega_pl0=0;

  /*********************
  variables in QV
  **********************/
  double interT=0, interL=0;
  double inter_imag_fxcT=0;
  double omega_tilde=0;
  double *a_3_L=(double *)malloc(6*sizeof(double));
  double *a_3_T=(double *)malloc(6*sizeof(double));
  double *b_3_L=(double *)malloc(6*sizeof(double));
  double *b_3_T=(double *)malloc(6*sizeof(double));
  double *gamma_3_L=(double *)malloc(6*sizeof(double));
  double *gamma_3_T=(double *)malloc(6*sizeof(double));
  double *omega_3_L=(double *)malloc(6*sizeof(double));
  double *omega_3_T=(double *)malloc(6*sizeof(double));
  
  /***************************************
  For integral in Kramers-Kronig relation
  ****************************************/ 
  double omega_min=0.0;
  double omega_max=0.0; 
  double ho=0.0;
  double omega_epsilon=0.0;
  
  for(ii=0;ii<=5;ii++){
     n[ii]=0; mu[ii]=0; f_T[ii]=0;
     k_F[ii]=0; epsilon_F[ii]=0; omega_pl[ii]=0;
     fxcL_0[ii]=0; fxcL_inf[ii]=0; fxcT_inf[ii]=0;
     beta[ii]=0; c0[ii]=0; c1[ii]=0;
     omega1[ii]=0; omega2[ii]=0; d0[ii]=0; d1[ii]=0;
     imag_fxcL[ii]=0; imag_fxcT[ii]=0;
     real_fxcL[ii]=0; real_fxcT[ii]=0;
     a_3_L[ii]=0; a_3_T[ii]=0;
     b_3_L[ii]=0; b_3_T[ii]=0;
     gamma_3_L[ii]=0; gamma_3_T[ii]=0;
     omega_3_L[ii]=0; omega_3_T[ii]=0; 
  }

  if(icase==1){
    printf("\n\t *** NCT0 ***\n");
    /**************************************************
    Calculate the values of f_xcT on r_s=0.5, 1,...,5
    needed in cubic spline interpolation for f_xcT or 
    other higher order interpolation
    **************************************************/
   
    mu[0]=sqrt(27.0)/sqrt(pow(0.5,9.0))/2.0/PI*0.0065;  // r_s=0.5
    mu[1]=sqrt(27.0)/sqrt(pow(1.0,9.0))/2.0/PI*0.0064;  // r_s=1
    mu[2]=sqrt(27.0)/sqrt(pow(2.0,9.0))/2.0/PI*0.0052;  // r_s=2
    mu[3]=sqrt(27.0)/sqrt(pow(3.0,9.0))/2.0/PI*0.0037;  // r_s=3
    mu[4]=sqrt(27.0)/sqrt(pow(4.0,9.0))/2.0/PI*0.0020;  // r_s=4
    mu[5]=sqrt(27.0)/sqrt(pow(5.0,9.0))/2.0/PI*0.0002;  // r_s=5

    n[0]=1.0/(4.0*PI/3.0*pow(0.5,3.0));  // r_s=0.5
    n[1]=1.0/(4.0*PI/3.0*pow(1.0,3.0));  // r_s=1
    n[2]=1.0/(4.0*PI/3.0*pow(2.0,3.0));  // r_s=2  
    n[3]=1.0/(4.0*PI/3.0*pow(3.0,3.0));  // r_s=3  
    n[4]=1.0/(4.0*PI/3.0*pow(4.0,3.0));  // r_s=4  
    n[5]=1.0/(4.0*PI/3.0*pow(5.0,3.0));  // r_s=5

    /**************************************
    nodes values at r_s=0.5, 1, 2, 3, 4, 5
    ***************************************/
    for(ii=0;ii<=5;ii++){
       f_T[ii]=mu[ii]/(n[ii]*n[ii]); 
       //printf("%f\n",f_T[ii]);
    }  
    //getchar();
    /******************
    calculate f_xcL,T
    *******************/
    for(ix=0;ix<Nx;ix++){
       for(iy=0;iy<Ny;iy++){
          for(iz=0;iz<Nz;iz++){  
             f_xcL[iz+Nz*(ix+Nx*iy)]=0;  

             r_s=cabs(pow(3.0/4.0/PI/rho_0[iz+Nz*(ix+Nx*iy)],1.0/3.0));  

             if((r_s)>=5){
               f_xcT[iz+Nz*(ix+Nx*iy)]=f_T[5];
             }                       
             else{ 
               /************************************
               4-th degree polynomial use f1,...,f5
               ************************************/
               f_xcT[iz+Nz*(ix+Nx*iy)]=(r_s-x2)*(r_s-x3)*(r_s-x4)*(r_s-x5)/(24.0)*f_T[1]
                                      +(r_s-x1)*(r_s-x3)*(r_s-x4)*(r_s-x5)/(-6.0)*f_T[2]
                                      +(r_s-x1)*(r_s-x2)*(r_s-x4)*(r_s-x5)/(4.0)*f_T[3]
                                      +(r_s-x1)*(r_s-x2)*(r_s-x3)*(r_s-x5)/(-6.0)*f_T[4]
                                      +(r_s-x1)*(r_s-x2)*(r_s-x3)*(r_s-x4)/(24.0)*f_T[5];     

  
             }  
          }
       }
    }
  }
  else if(icase==2){
    printf("\n\t *** NCT ***\n");
    /******
    NCT
    *******/
    /**************************************************
    Calculate the values of f_xcT,L on r_s=0.5, 1,...,5
    needed in cubic spline interpolation for f_xcT,L or 
    other higher order interpolation
    **************************************************/ 
    n[0]=1.0/(4.0*PI/3.0*pow(0.5,3.0));  // r_s=0.5
    n[1]=1.0/(4.0*PI/3.0*pow(1.0,3.0));  // r_s=1
    n[2]=1.0/(4.0*PI/3.0*pow(2.0,3.0));  // r_s=2  
    n[3]=1.0/(4.0*PI/3.0*pow(3.0,3.0));  // r_s=3  
    n[4]=1.0/(4.0*PI/3.0*pow(4.0,3.0));  // r_s=4  
    n[5]=1.0/(4.0*PI/3.0*pow(5.0,3.0));  // r_s=5

    for(ii=0;ii<=5;ii++){
       omega_pl[ii]=sqrt(4.0*PI*n[ii]);       // plasmon frequency
       k_F[ii]=pow(3.0*PI*PI*n[ii],1.0/3.0);  // Fermi momentum
       epsilon_F[ii]=k_F[ii]*k_F[ii]/2.0/omega_pl[ii];     // Fermi energy
       //printf("%f,%f,%f\n",omega_pl[ii], k_F[ii], epsilon_F[ii]);
    }
    //getchar();

    /******************************************************
    parameters in analytical interpolation for imag_fxcL,T
    *******************************************************/
    fxcL_0[0]=-0.04246;
    fxcL_0[1]=-0.0611;
    fxcL_0[2]=-0.0891;
    fxcL_0[3]=-0.1119;
    fxcL_0[4]=-0.1320;
    fxcL_0[5]=-0.1503;

    fxcL_inf[0]=-0.01794;
    fxcL_inf[1]=-0.0216;
    fxcL_inf[2]=-0.0252;
    fxcL_inf[3]=-0.0280;
    fxcL_inf[4]=-0.0308;
    fxcL_inf[5]=-0.0338;

    fxcT_inf[0]=0.0177;
    fxcT_inf[1]=0.0284;
    fxcT_inf[2]=0.0457;
    fxcT_inf[3]=0.0600;
    fxcT_inf[4]=0.0724;
    fxcT_inf[5]=0.0835;

    beta[0]=1.87;
    beta[1]=1.48;
    beta[2]=1.22;
    beta[3]=1.1;
    beta[4]=1.02;
    beta[5]=0.955;

    c0[0]=0.175/100.0;
    c0[1]=0.421/100.0;
    c0[2]=0.895/100.0;
    c0[3]=1.29/100.0;
    c0[4]=1.65/100.0;
    c0[5]=1.94/100.0;

    c1[0]=0.694/100.0;
    c1[1]=1.76/100.0;
    c1[2]=3.87/100.0;
    c1[3]=6.09/100.0;
    c1[4]=7.87/100.0;
    c1[5]=9.82/100.0;

    omega1[0]=1.75;
    omega1[1]=0.982;
    omega1[2]=0.347;
    omega1[3]=0.143;
    omega1[4]=-0.143;
    omega1[5]=-0.27;

    omega2[0]=-3.59;
    omega2[1]=-1.45;
    omega2[2]=0.181;
    omega2[3]=0.693;
    omega2[4]=1.33;
    omega2[5]=1.61;

    d0[0]=0.173;
    d0[1]=0.291;
    d0[2]=0.49;
    d0[3]=0.664;
    d0[4]=0.824;
    d0[5]=0.974;

    d1[0]=5.72/100.0;
    d1[1]=9.38/100.0;
    d1[2]=13.2/100.0;
    d1[3]=16.7/100.0;
    d1[4]=17.0/100.0;
    d1[5]=18.3/100.0;

    /************************************* 
    first, find f_xcL,T at r_s=0.5,1,...,5
    *************************************/
    for(ii=0;ii<=5;ii++){
       /****************
       Imaginary part 
       *****************/
       omega0=omega/omega_pl[ii];  
       if(omega0>=0){
         g_x=(beta[ii]+0.5*omega0/2.0/epsilon_F[ii])/(1.0+omega0/2.0/epsilon_F[ii]);
         if(omega0<=2){
           imag_fxcL[ii]=-g_x*(c0[ii]*omega0+c1[ii]*(omega0-1.0)/(exp(7.0/omega0-5.0)+1.0));
         }
         else{
           imag_fxcL[ii]=-g_x*(d0[ii]*sqrt(omega0-2.0)+d1[ii])/(omega0*(omega0-omega1[ii]*sqrt(omega0)-omega2[ii]));
         }
         imag_fxcT[ii]=0.72*imag_fxcL[ii]; 
       }
       else{ 
         g_x=(beta[ii]+0.5*fabs(omega0)/2.0/epsilon_F[ii])/(1.0+fabs(omega0)/2.0/epsilon_F[ii]);
         if(fabs(omega0)<=2){
           imag_fxcL[ii]=-g_x*(c0[ii]*fabs(omega0)+c1[ii]*(fabs(omega0)-1.0)/(exp(7.0/fabs(omega0)-5.0)+1.0));
         }
         else{
           imag_fxcL[ii]=-g_x*(d0[ii]*sqrt(fabs(omega0)-2.0)+d1[ii])/(fabs(omega0)*(fabs(omega0)-omega1[ii]*sqrt(fabs(omega0))-omega2[ii]));
         }
         imag_fxcL[ii]=-imag_fxcL[ii];
         imag_fxcT[ii]=0.72*imag_fxcL[ii]; 
       }
       
       /******************
       Real part
       ******************/
       omega_min=omega0-100.0; 
       omega_max=omega0+100.0;
       omega_epsilon=0.0001; ho=0.001; 
       for(omegaP=omega0-omega_epsilon; omegaP>=omega_min; omegaP=omegaP-ho){
          if(omegaP>=0){ 
            g_x=(beta[ii]+0.5*omegaP/2.0/epsilon_F[ii])/(1.0+omegaP/2.0/epsilon_F[ii]);
            if(omegaP<=2){
              inter_imag_fxcL=-g_x*(c0[ii]*omegaP+c1[ii]*(omegaP-1.0)/(exp(7.0/omegaP-5.0)+1.0));
            }
            else{
              inter_imag_fxcL=-g_x*(d0[ii]*sqrt(omegaP-2.0)+d1[ii])/(omegaP*(omegaP-omega1[ii]*sqrt(omegaP)-omega2[ii]));
            }
            inter=inter+inter_imag_fxcL/(omegaP-omega0)*ho;
          }
          else{ 
            g_x=(beta[ii]+0.5*fabs(omegaP)/2.0/epsilon_F[ii])/(1.0+fabs(omegaP)/2.0/epsilon_F[ii]);
            if(fabs(omegaP)<=2){
              inter_imag_fxcL=-g_x*(c0[ii]*fabs(omegaP)+c1[ii]*(fabs(omegaP)-1.0)/(exp(7.0/fabs(omegaP)-5.0)+1.0));
            }
            else{
              inter_imag_fxcL=-g_x*(d0[ii]*sqrt(fabs(omegaP)-2.0)+d1[ii])/(fabs(omegaP)*(fabs(omegaP)-omega1[ii]*sqrt(fabs(omegaP))-omega2[ii]));
            }
            inter_imag_fxcL=-inter_imag_fxcL;
            inter=inter+inter_imag_fxcL/(omegaP-omega0)*ho;
          } 
       }

       for(omegaP=omega0+omega_epsilon; omegaP<=omega_max; omegaP=omegaP+ho){
          if(omegaP>=0){ 
            g_x=(beta[ii]+0.5*omegaP/2.0/epsilon_F[ii])/(1.0+omegaP/2.0/epsilon_F[ii]);
            if(omegaP<=2){
              inter_imag_fxcL=-g_x*(c0[ii]*omegaP+c1[ii]*(omegaP-1.0)/(exp(7.0/omegaP-5.0)+1.0));
            }
            else{
              inter_imag_fxcL=-g_x*(d0[ii]*sqrt(omegaP-2.0)+d1[ii])/(omegaP*(omegaP-omega1[ii]*sqrt(omegaP)-omega2[ii]));
            }
            inter=inter+inter_imag_fxcL/(omegaP-omega0)*ho;
          }
          else{ 
            g_x=(beta[ii]+0.5*fabs(omegaP)/2.0/epsilon_F[ii])/(1.0+fabs(omegaP)/2.0/epsilon_F[ii]);
            if(fabs(omegaP)<=2){
              inter_imag_fxcL=-g_x*(c0[ii]*fabs(omegaP)+c1[ii]*(fabs(omegaP)-1.0)/(exp(7.0/fabs(omegaP)-5.0)+1.0));
            }
            else{
              inter_imag_fxcL=-g_x*(d0[ii]*sqrt(fabs(omegaP)-2.0)+d1[ii])/(fabs(omegaP)*(fabs(omegaP)-omega1[ii]*sqrt(fabs(omegaP))-omega2[ii]));
            }
            inter_imag_fxcL=-inter_imag_fxcL;
            inter=inter+inter_imag_fxcL/(omegaP-omega0)*ho;
          } 
       } 
       inter=inter/PI;    
       real_fxcL[ii]=fxcL_inf[ii]+inter;
       real_fxcT[ii]=fxcT_inf[ii]+0.72*inter;

       real_fxcL[ii]=real_fxcL[ii]*2.0*omega_pl[ii]/n[ii]; real_fxcT[ii]=real_fxcT[ii]*2.0*omega_pl[ii]/n[ii];  // in a.u.
       imag_fxcL[ii]=imag_fxcL[ii]*2.0*omega_pl[ii]/n[ii]; imag_fxcT[ii]=imag_fxcT[ii]*2.0*omega_pl[ii]/n[ii];  // in a.u.
    }
       
 
    /*****************************************************************************
    second, calculate f_xcL,T using interpolation with values at r_s=0.5, 1,...,5
    *****************************************************************************/
    for(ix=0;ix<Nx;ix++){
       for(iy=0;iy<Ny;iy++){
          for(iz=0;iz<Nz;iz++){  
             r_s=cabs(pow(3.0/4.0/PI/rho_0[iz+Nz*(ix+Nx*iy)],1.0/3.0)); 
             //omega_pl0=sqrt(4.0*PI*rho_0[iz+Nz*(ix+Nx*iy)]);

             if((r_s)>=5){
               f_xcT[iz+Nz*(ix+Nx*iy)]=real_fxcT[5]+_Complex_I*imag_fxcT[5];
               f_xcL[iz+Nz*(ix+Nx*iy)]=real_fxcL[5]+_Complex_I*imag_fxcL[5];
             }
             else{ 
               /************************************
               4-th degree polynomial use f1,...,f5
               ************************************/
               f_xcT[iz+Nz*(ix+Nx*iy)]=(r_s-x2)*(r_s-x3)*(r_s-x4)*(r_s-x5)/(24.0)*(real_fxcT[1]+_Complex_I*imag_fxcT[1])
                                      +(r_s-x1)*(r_s-x3)*(r_s-x4)*(r_s-x5)/(-6.0)*(real_fxcT[2]+_Complex_I*imag_fxcT[2])
                                      +(r_s-x1)*(r_s-x2)*(r_s-x4)*(r_s-x5)/(4.0)* (real_fxcT[3]+_Complex_I*imag_fxcT[3])
                                      +(r_s-x1)*(r_s-x2)*(r_s-x3)*(r_s-x5)/(-6.0)*(real_fxcT[4]+_Complex_I*imag_fxcT[4])
                                      +(r_s-x1)*(r_s-x2)*(r_s-x3)*(r_s-x4)/(24.0)*(real_fxcT[5]+_Complex_I*imag_fxcT[5]);  

               
               f_xcL[iz+Nz*(ix+Nx*iy)]=(r_s-x2)*(r_s-x3)*(r_s-x4)*(r_s-x5)/(24.0)*(real_fxcL[1]+_Complex_I*imag_fxcL[1])
                                      +(r_s-x1)*(r_s-x3)*(r_s-x4)*(r_s-x5)/(-6.0)*(real_fxcL[2]+_Complex_I*imag_fxcL[2])
                                      +(r_s-x1)*(r_s-x2)*(r_s-x4)*(r_s-x5)/(4.0)* (real_fxcL[3]+_Complex_I*imag_fxcL[3])
                                      +(r_s-x1)*(r_s-x2)*(r_s-x3)*(r_s-x5)/(-6.0)*(real_fxcL[4]+_Complex_I*imag_fxcL[4])
                                      +(r_s-x1)*(r_s-x2)*(r_s-x3)*(r_s-x4)/(24.0)*(real_fxcL[5]+_Complex_I*imag_fxcL[5]);      

  
             }   
          }
       }
    } 
  }
  else if(icase==3){
    printf("\t *** QV ***\n");
    /*****
    QV
    *****/ 
    n[0]=1.0/(4.0*PI/3.0*pow(0.5,3.0));  // r_s=0.5
    n[1]=1.0/(4.0*PI/3.0*pow(1.0,3.0));  // r_s=1
    n[2]=1.0/(4.0*PI/3.0*pow(2.0,3.0));  // r_s=2  
    n[3]=1.0/(4.0*PI/3.0*pow(3.0,3.0));  // r_s=3  
    n[4]=1.0/(4.0*PI/3.0*pow(4.0,3.0));  // r_s=4  
    n[5]=1.0/(4.0*PI/3.0*pow(5.0,3.0));  // r_s=5 
    /**********************
    parameters
    ***********************/
    /* note that fxcL,T_inf is in unit 2w_pl/n */ 
    fxcL_inf[0]=-0.01794;
    fxcL_inf[1]=-0.0216;
    fxcL_inf[2]=-0.0252;
    fxcL_inf[3]=-0.0280;
    fxcL_inf[4]=-0.0308;
    fxcL_inf[5]=-0.0338;

    fxcT_inf[0]=0.0177;
    fxcT_inf[1]=0.0284;
    fxcT_inf[2]=0.0457;
    fxcT_inf[3]=0.0600;
    fxcT_inf[4]=0.0724;
    fxcT_inf[5]=0.0835; 

    a_3_L[0]=0; 
    a_3_L[1]=0.5026/100.0; 
    a_3_L[2]=0.8473/100.0; 
    a_3_L[3]=1.092/100.0; 
    a_3_L[4]=1.278/100.0; 
    a_3_L[5]=1.426/100.0; 
    
    b_3_L[0]=0; 
    b_3_L[1]=0.1555;   
    b_3_L[2]=0.1558;     
    b_3_L[3]=0.1496;    
    b_3_L[4]=0.1428;    
    b_3_L[5]=0.1363;
    
    gamma_3_L[0]=0; 
    gamma_3_L[1]=1.656;  
    gamma_3_L[2]=1.368;  
    gamma_3_L[3]=1.215;   
    gamma_3_L[4]=1.112;   
    gamma_3_L[5]=1.033;
    
    omega_3_L[0]=0; 
    omega_3_L[1]=-1.484; 
    omega_3_L[2]=-1.052; 
    omega_3_L[3]=-0.8227; 
    omega_3_L[4]=-0.6683; 
    omega_3_L[5]=-0.5498;
    
    a_3_T[0]=0; 
    a_3_T[1]=0.3769/100.0; 
    a_3_T[2]=0.6355/100.0; 
    a_3_T[3]=0.8191/100.0; 
    a_3_T[4]=0.9587/100.0; 
    a_3_L[5]=1.069/100.0; 
    
    b_3_T[0]=0; 
    b_3_T[1]=0.1651;   
    b_3_T[2]=0.1654;     
    b_3_T[3]=0.1589;    
    b_3_T[4]=0.1516;    
    b_3_T[5]=0.1448;
    
    gamma_3_T[0]=0; 
    gamma_3_T[1]=1.821;  
    gamma_3_T[2]=1.533;  
    gamma_3_T[3]=1.380;   
    gamma_3_T[4]=1.277;   
    gamma_3_T[5]=1.198;
 
    omega_3_T[0]=0; 
    omega_3_T[1]=-1.732; 
    omega_3_T[2]=-1.300; 
    omega_3_T[3]=-1.070;  
    omega_3_T[4]=-0.9158; 
    omega_3_T[5]=-0.7973;


    /*******************************
    f_xcL,T at r_s=1,..,5
    ********************************/ 
    for(ii=0;ii<5;ii++){
       /***************
       Imaginary part
       ***************/
       omega_pl[ii]=sqrt(4.0*PI*n[ii]);       // in a.u.
       omega_tilde=omega/(2.0*omega_pl[ii]);  // in a.u.
       imag_fxcL[ii]=-2.0*omega_pl[ii]/n[ii]*omega_tilde*(a_3_L[ii]/pow(1.0+b_3_L[ii]*omega_tilde*omega_tilde,5.0/4.0)+omega_tilde*omega_tilde*cexp(-pow(fabs(omega_tilde)-omega_3_L[ii],2.0)/gamma_3_L[ii]));
       imag_fxcT[ii]=-2.0*omega_pl[ii]/n[ii]*omega_tilde*(a_3_T[ii]/pow(1.0+b_3_T[ii]*omega_tilde*omega_tilde,5.0/4.0)+omega_tilde*omega_tilde*cexp(-pow(fabs(omega_tilde)-omega_3_T[ii],2.0)/gamma_3_T[ii]));

       /******************
       Real part
       ******************/
       omega_min=omega-1000.0; 
       omega_max=omega+1000.0;
       omega_epsilon=0.0001; 
       ho=0.01; 
       interL=0; interT=0;
       for(omegaP=omega-omega_epsilon; omegaP>=omega_min; omegaP=omegaP-ho){  
          omega_tilde=omegaP/(2.0*omega_pl[ii]);
          inter_imag_fxcL=-2.0*omega_pl[ii]/n[ii]*omega_tilde*(a_3_L[ii]/pow(1.0+b_3_L[ii]*omega_tilde*omega_tilde,5.0/4.0)+omega_tilde*omega_tilde*cexp(-pow(fabs(omega_tilde)-omega_3_L[ii],2.0)/gamma_3_L[ii]));

          inter_imag_fxcT=-2.0*omega_pl[ii]/n[ii]*omega_tilde*(a_3_T[ii]/pow(1.0+b_3_T[ii]*omega_tilde*omega_tilde,5.0/4.0)+omega_tilde*omega_tilde*cexp(-pow(fabs(omega_tilde)-omega_3_T[ii],2.0)/gamma_3_T[ii])); 
          
          interL=interL+inter_imag_fxcL/(omegaP-omega)*ho;
          interT=interT+inter_imag_fxcT/(omegaP-omega)*ho;
       }
       for(omegaP=omega+omega_epsilon; omegaP<=omega_max; omegaP=omegaP+ho){  
          omega_tilde=omegaP/(2.0*omega_pl[ii]);
          inter_imag_fxcL=-2.0*omega_pl[ii]/n[ii]*omega_tilde*(a_3_L[ii]/pow(1.0+b_3_L[ii]*omega_tilde*omega_tilde,5.0/4.0)+omega_tilde*omega_tilde*cexp(-pow(fabs(omega_tilde)-omega_3_L[ii],2.0)/gamma_3_L[ii]));

          inter_imag_fxcT=-2.0*omega_pl[ii]/n[ii]*omega_tilde*(a_3_T[ii]/pow(1.0+b_3_T[ii]*omega_tilde*omega_tilde,5.0/4.0)+omega_tilde*omega_tilde*cexp(-pow(fabs(omega_tilde)-omega_3_T[ii],2.0)/gamma_3_T[ii])); 
          
          interL=interL+inter_imag_fxcL/(omegaP-omega)*ho;
          interT=interT+inter_imag_fxcT/(omegaP-omega)*ho;
       }
       interL=interL/PI;
       interT=interT/PI;       
       real_fxcL[ii]=fxcL_inf[ii]*(2.0*omega_pl[ii]/n[ii])+interL;  // in a.u.
       real_fxcT[ii]=fxcT_inf[ii]*(2.0*omega_pl[ii]/n[ii])+interT;  // in a.u. 
    } 

    /********************************
    Interpolation to all r_s
    *********************************/
    for(ix=0;ix<Nx;ix++){
       for(iy=0;iy<Ny;iy++){
          for(iz=0;iz<Nz;iz++){  
             r_s=cabs(pow(3.0/4.0/PI/rho_0[iz+Nz*(ix+Nx*iy)],1.0/3.0));  

             if((r_s)>=5){
               f_xcT[iz+Nz*(ix+Nx*iy)]=real_fxcT[5]+_Complex_I*imag_fxcT[5];
               f_xcL[iz+Nz*(ix+Nx*iy)]=real_fxcL[5]+_Complex_I*imag_fxcL[5];
             }
             else{ 
               /************************************
               4-th degree polynomial use f1,...,f5
               ************************************/
               f_xcT[iz+Nz*(ix+Nx*iy)]=(r_s-x2)*(r_s-x3)*(r_s-x4)*(r_s-x5)/(24.0)*(real_fxcT[1]+_Complex_I*imag_fxcT[1])
                                      +(r_s-x1)*(r_s-x3)*(r_s-x4)*(r_s-x5)/(-6.0)*(real_fxcT[2]+_Complex_I*imag_fxcT[2])
                                      +(r_s-x1)*(r_s-x2)*(r_s-x4)*(r_s-x5)/(4.0)* (real_fxcT[3]+_Complex_I*imag_fxcT[3])
                                      +(r_s-x1)*(r_s-x2)*(r_s-x3)*(r_s-x5)/(-6.0)*(real_fxcT[4]+_Complex_I*imag_fxcT[4])
                                      +(r_s-x1)*(r_s-x2)*(r_s-x3)*(r_s-x4)/(24.0)*(real_fxcT[5]+_Complex_I*imag_fxcT[5]);  

               
               f_xcL[iz+Nz*(ix+Nx*iy)]=(r_s-x2)*(r_s-x3)*(r_s-x4)*(r_s-x5)/(24.0)*(real_fxcL[1]+_Complex_I*imag_fxcL[1])
                                      +(r_s-x1)*(r_s-x3)*(r_s-x4)*(r_s-x5)/(-6.0)*(real_fxcL[2]+_Complex_I*imag_fxcL[2])
                                      +(r_s-x1)*(r_s-x2)*(r_s-x4)*(r_s-x5)/(4.0)* (real_fxcL[3]+_Complex_I*imag_fxcL[3])
                                      +(r_s-x1)*(r_s-x2)*(r_s-x3)*(r_s-x5)/(-6.0)*(real_fxcL[4]+_Complex_I*imag_fxcL[4])
                                      +(r_s-x1)*(r_s-x2)*(r_s-x3)*(r_s-x4)/(24.0)*(real_fxcL[5]+_Complex_I*imag_fxcL[5]);      

  
             }   
          }
       }
    }

  }
  else{
    printf("Wrong choice of f_xcL,T \n"); getchar();
  }

  free(mu); free(n); free(f_T);

  free(k_F); free(epsilon_F); free(omega_pl);

  free(fxcL_0); free(fxcL_inf); free(fxcT_inf);
  free(beta); free(c0); free(c1); 
  free(omega1); free(omega2); 
  free(d0); free(d1);

  free(imag_fxcL); free(imag_fxcT);
  free(real_fxcL); free(real_fxcT); 

  free(a_3_L); free(a_3_T);
  free(b_3_L); free(b_3_T);
  free(gamma_3_L); free(gamma_3_T);
  free(omega_3_L); free(omega_3_T);
}

/**
 * 3.
 * 
 */

void f_xc_ALDA(double complex *f_xc, 
	       double *rho_0, double *QD_epsilon, 
	       int Nx, int Ny, int Nz)
{
  int ix=0, iy=0, iz=0;

  /*************************
  Parameter in VWN
  **************************/
  double A=0, x_0=0, b=0, c=0, x=0, X=0, Q=0, X_0=0, epsilon_x=0, epsilon_c=0, de_x_dr=0, dx_dr=0, dX=0, de_c_dr=0, ddX=0, dde_x_ddr=0, dde_c_ddr=0;
  double de_c_dx=0, ddx_ddr=0, dde_c_ddx=0;;
  //double rho_02=0;

  /***************************
  Parameter in Hedin-Lundqvist
  ****************************/ 
  int id;
  double rs, beta, exc, drdn, ddrddn, dbetadr, ddbetaddr, dbetadn, ddbetaddn, dexcdn, ddexcddn;

  /***************************
  Parameter in Perdew-Zunger
  ****************************/ 
  double dexdrho, drdrho, decdr, decdrho, dexcdrho, ddexddrho, ddrddrho, ddecddr, ddecddrho;   
  double gamma0=-0.1423, beta1=1.0529, beta2=0.3334, A0=0.0311, B0=-0.0480, C0=0.0020, D0=-0.0116;

  int icase=1;
  /*************************************
  icase=1, LDA  VWN
  icase=2, LDA  Hedin-Lundqvist
  icase=3, LDA  Perdew-Zunger
  **************************************/
  if(icase==1){
    printf("\t ALDA with VMN ******\n");
  }
  else if(icase==2){
    printf("\t ALDA with Hdein-Lundqvist ******\n");
  }
  else if(icase==3){
    printf("\t ALDA with Perdew-Zunger ******\n");
  }
  else{
    printf("\t ALDA with wrong choice, use default VMN******\n");
    icase=1;
  }

  if(icase==1){
    A=0.0310907; x_0=-0.10498; b=3.72744; c=12.9352; Q=sqrt(4.0*c-b*b); X_0=x_0*x_0+b*x_0+c; ddX=2.0;
    for(ix=0;ix<Nx;ix++){
       for(iy=0;iy<Ny;iy++){
          for(iz=0;iz<Nz;iz++){
             id=iz+Nz*(ix+Nx*iy);
             //rho_02=rho_0[id]*rho_0[id]; 
             
             x=pow(3.0/4.0/PI/rho_0[id], 1.0/6.0);

             X=x*x+b*x+c;
             
             de_x_dr=-1.0/4.0/PI*pow(3.0*PI*PI,1.0/3.0)/pow(rho_0[id], 2.0/3.0);

             dx_dr=-1.0/6.0*pow(3.0/4.0/PI,1.0/6.0)/pow(rho_0[id],7.0/6.0);

             dX=2.0*x+b;

             de_c_dx=(A*((2.0*X-x*dX)/(x*X)-4.0*b/((2.0*x+b)*(2.0*x+b)+Q*Q)-b*x_0/X_0 *((2.0*X-(x-x_0)*dX)/((x-x_0)*X)-4.0*(b+2.0*x_0)/((2.0*x+b)*(2.0*x+b)+Q*Q))));
             de_c_dr=de_c_dx*dx_dr;
  
             ddx_ddr=7.0/36.0*pow(3.0/4.0/PI,1.0/6.0)/pow(rho_0[id],13.0/6.0);

             dde_x_ddr=1.0/6.0/PI*pow(3.0*PI*PI,1.0/3.0)/pow(rho_0[id],5.0/3.0); 

             dde_c_ddx=(A*((dX-x*ddX)/(x*X)-(2.0*X-x*dX)*(X+x*dX)/(x*x*X*X)+16.0*b*(2.0*x+b)/(((2.0*x+b)*(2.0*x+b)+Q*Q)*((2.0*x+b)*(2.0*x+b)+Q*Q))-b*x_0/X_0*((dX-(x-x_0)*ddX)/((x-x_0)*X)-(2.0*X-(x-x_0)*dX)*(X+(x-x_0)*dX)/((x-x_0)*(x-x_0)*X*X)+16.0*(2.0*x_0+b)*(2.0*x+b)/(((2.0*x+b)*(2.0*x+b)+Q*Q)*((2.0*x+b)*(2.0*x+b)+Q*Q)))));

             dde_c_ddr=dde_c_ddx*dx_dr*dx_dr+de_c_dx*ddx_ddr;

             f_xc[id]=2.0*(de_x_dr+de_c_dr)+ rho_0[id]*(dde_x_ddr+dde_c_ddr);  f_xc[id]=f_xc[id]/QD_epsilon[id];
          }
       }
    }
  }
  else if(icase==2){
    for(ix=0;ix<Nx;ix++){
       for(iy=0;iy<Ny;iy++){
          for(iz=0;iz<Nz;iz++){
             id=iz+Nz*(ix+Nx*iy);
             rs=pow(3.0/4.0/PI/rho_0[id], 1.0/3.0);
             beta=1.0+0.0368*rs*log(1.0+21.0/rs);
             exc=-2.0*beta*pow(3.0*rho_0[id]/8.0/PI, 1.0/3.0);

             drdn=1.0/3.0*pow(3.0/4.0/PI,1.0/3.0)*pow(rho_0[id], -2.0/3.0);
             ddrddn=-2.0/9.0*pow(rho_0[id],-5.0/3.0)*pow(3.0/4.0/PI,1.0/3.0);
             dbetadr=0.0368*log(1.0+21.0/rs)-0.0368*21.0/(rs+21.0);
             ddbetaddr=-0.0368*21.0*21.0/rs/pow(rs+21.0,2.0);
             dbetadn=dbetadr*drdn;
             ddbetaddn=ddbetaddr*pow(drdn,2.0)+dbetadr*ddrddn;

             dexcdn=-2.0*pow(3.0/8.0/PI,1.0/3.0)*(dbetadn*pow(rho_0[id],1.0/3.0)+beta/3.0*pow(rho_0[id],-2.0/3.0));
             ddexcddn=-2.0*pow(3.0/8.0/PI,1.0/3.0)*(ddbetaddn*pow(rho_0[id],1.0/3.0)+2.0*dbetadn/3.0*pow(rho_0[id],-2.0/3.0)-beta*2.0/9.0*pow(rho_0[id],-5.0/3.0));

             f_xc[id]=ddexcddn*rho_0[id]+2.0*dexcdn;  f_xc[id]=f_xc[id]/QD_epsilon[id];
          }
       }
    }
  }
  else if(icase==3){
    for(ix=0;ix<Nx;ix++){
       for(iy=0;iy<Ny;iy++){
          for(iz=0;iz<Nz;iz++){
             id=iz+Nz*(ix+Nx*iy);
             rs=pow(3.0/4.0/PI/rho_0[id], 1.0/3.0);
             dexdrho=pow(3.0/PI,1.0/3.0)/4.0*pow(rho_0[id],-2.0/3.0);
             drdrho=-pow(3.0/4.0/PI,1.0/3.0)/3.0*pow(rho_0[id],-4.0/3.0);
             if(rs<1){
               decdr=A0/rs+C0*(log(rs)+1)+D0;
             }
             else{
               decdr=-gamma0*(beta1/2.0/sqrt(rs)+beta2)/pow(1+beta1*sqrt(rs)+beta2*rs,2.0);
             } 

             decdrho=decdr*drdrho;

             dexcdrho=dexdrho+decdrho;
             
             ddexddrho=-pow(3.0/PI,1.0/3.0)/6.0*pow(rho_0[id],-5.0/3.0);

             ddrddrho=4.0/9.0*pow(3.0/4.0/PI,1.0/3.0)*pow(rho_0[id],-7.0/3.0);

             if(rs<1){
               ddecddr=-A0/rs/rs+C0/rs;
             }
             else{
               ddecddr=-gamma0*(-beta1/4.0*pow(rho_0[id],-3.0/2.0)*(1+beta1*sqrt(rs)+beta2*rs)-2.0*pow(beta1/2.0/sqrt(rs)+beta2,2.0))/pow(1+beta1*sqrt(rs)+beta2*rs,3.0);
             }

             ddecddrho=ddecddr*drdrho*drdrho+decdr*ddrddrho; 

             f_xc[id]=2*dexcdrho+rho_0[id]*(ddexddrho+ddecddrho);    f_xc[id]=f_xc[id]/QD_epsilon[id];
          }
       }
    }
  }
  else{
    printf("Wrong reference in eta_xc_xi_xc\n"); getchar();
  }

}

/**
 * 
 * 
 */

void eta_xc_xi_xc(double complex *eta_xc, double complex *xi_xc, 
		  double complex *fxcL0, double complex *f_xcL, double complex *f_xcT, 
		  double *rho_0, double complex omega_input, int Nx, int Ny, int Nz)
{
  int ix=0, iy=0, iz=0, idx; 
  double complex rho_02=0;
  //double omega=creal(omega_input);
  double complex omega = omega_input;
  for(iy=0;iy<Ny;iy++){
    for(ix=0;ix<Nx;ix++){
      for(iz=0;iz<Nz;iz++){
	idx = iz + Nz*(ix + Nx*iy);
	rho_02 = rho_0[idx]*rho_0[idx];
	eta_xc[idx] = -rho_02/omega/_Complex_I*f_xcT[idx]; 
	xi_xc[idx] = -rho_02/omega/_Complex_I*
	  (f_xcL[idx] - 4.0/3.0*f_xcT[idx] - fxcL0[idx]); 
      }
    }
  }  

}
